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ABSTRACT 



We use cosmological simulations of reionization to predict the effect of large H II regions 
around individual high-redshift quasars on the possible signal from the redshifted 21cm line of 
neutral hydrogen in the pre-reionization era. We show that these H II regions appear as "spectral 
dips" in frequency space with equivalent widths in excess of 3 mK MHz or depths in excess of 
about 1.5 mK, and that they are by far the most prominent cosmological signals in the redshifted 
H I distribution. 

These spectral dips are expected to be present in almost every line of sight. If the spectral dips 
of a large enough sample of H II regions are well resolved in frequency space, the distribution of 
line depth and equivalent width in frequency with a known observing beamsize can be used to infer 
the H II region size distribution and the mean difference in neutral hydrogen density between the 
H II regions ( which may contain self-shielded neutral gas clumps) and the surrounding medium, 
providing a powerful test for models of reionization. 

Subject headings: cosmology: theory - cosmology: large-scale structure of universe - diffuse radiation - 
galaxies: formation - galaxies: intergalactic medium - radio lines: general 



1. Introduction 

Theoretical interest in the redshifted 21 cm line 
signal from intergalactic hydrogen before and dur- 
ing the reionization era has surged in recent years. 
Fluctuations in the 21 cm signal can be generated 
by variations in the density, spin temperature, or 
ionized fraction of the intergalactic gas. These 
could arise from primordial density fluctuations, 
from spatial variations of the radiation intensity 
near the wavelength of the Ly-a hydrogen line 
(which affects the spin temperature), from other 
radiation that can heat the gas, or from H II re- 
gions created during the process of reionization 



by discrete sources. The latter source of 21cm 
fluctuations was discussed in the pioneering work 
of Scott & Rees (1990), Madau, Meiksin, & Rees 
(1997), and Tozzi ct al. (2000), but has received 
only limited attention recently (Wyithe & Loeb 
2004). 

However, ionized (H II) regions around high- 
redshift quasars can easily reach sizes comparable 
to the angular beams of future radio telescopes 
such as LOFAR and SKA. Thus, they could be 
the most prominent features in the angular and 
frequency distribution of the redshifted 21 cm sig- 
nal from the early universe, and the first features 



1 



to be measured observationally. 

In this paper we estimate the effects of these 
large features on the possible cosmological sig- 
nal by using numerical simulations of reioniza- 
tion. This becomes feasible now as simulations 
of reionization reach the level of sophistication 
and accuracy sufficient to model quantitatively 
not only the small-scale details of reionization 
(Gnedin 2004), but also very large regions (up to 
a horizon size: Kohler, Gnedin, & Miralda-Escudc 
2004). The latter approach is used here, because 
quasars are relatively rare at high redshift, and 
sufficiently large computational volumes (several 
hundred Mpc in size) must be used to include a 
representative sample of quasars. 

We concentrate solely on the redshifted 21cm 
signal in the frequency domain, because, as has 
been shown before, foreground contamination is 
a severe limiting factor in observing the cosmo- 
logical features in the angular domain on the sky 
(Di Matteo et al. 2002; Oh & Mack 2004; Gnedin 
& Shaver 2004). Several recent papers have dis- 
cussed elegant and sophisticated approaches to 
circumventing this contamination (Zaldarriaga, 
Furlanctto & Hernquist 2004; He et al. 2004; San- 
tos et al. 2004, Locb & Zaldarriaga 2004), but 
these approaches would require bigger telescopes 
than currently planned, so it is most likely that 
the first detections will be made in the frequency 
domain and attempts to de-contaminate the an- 
gular signal will require follow-on observations. 

2. Simulation 

In order to model the spectral signatures of H II 
regions in the early universe, we need to simulate 
a large enough region of space that includes bright 
quasars during the reionization era. Such quasars 
are sparse at high redshift, since in the hierarchical 
structure formation paradigm they only form in 
extreme over-densities at early times. 

The simulation code we use is a modified ver- 
sion of the "Softened-Lagrangian Hydrodynamics" 
code described in Gnedin (2004), which follows the 
evolution of dark matter and gas. The code in- 
cludes a model for star formation in regions of 
dense gas that can cool, and for the emission of 
ionizing radiation from stars. Radiative transfer 
is fully included using the method described in 
Gnedin & Abel (2001), and the ionization state of 



the gas is followed in every cell as determined by 
the local radiation intensity. 

In this specific project we use a multi-resolution 
approach that utilizes a small-scale simulation as a 
"sub-cell" model for larger-scale simulations. The 
"sub-cell" model is implemented by using the so- 
called "clumping factors", which account for the 
structure on scales too small to be resolved in a 
large-scale simulation. The complete description 
of our approach and the appropriate tests are pre- 
sented in a separate paper (Kohler et al. 2004). 

Using the clumping factor formalism allows us 
to increase the size of the simulation box to sev- 
eral hundred comoving Mpc, while simultaneously 
accounting for the ionizing radiation transfer and 
recombination rate in a gaseous medium that is 
highly inhomogeneous on small scales. 

In our simulation of reionization, we add in- 
dividual quasars as discrete, luminous sources of 
ionizing radiation in addition to the emission asso- 
ciated with star-formation. The luminosity distri- 
bution of quasars added to the simulation is cho- 
sen from the model of Schirber & Bullock (2003). 
This model is consistent with all known observa- 
tional data on the evolution of the quasar luminos- 
ity function, including recent GOODS data (Cris- 
tiani et al. 2004). The impact of quasars on reion- 
ization depends of course on the luminosity func- 
tion at very high redshift, z — 6, before reioniza- 
tion ended, so we have to extrapolate the Schirber 
& Bullock model to redhifts higher than observed 
quasar redshifts, which introduces a substantial 
uncertainty. Our goal, therefore, is to discuss a 
range of possibilities, as precise predictions are not 
possible. The extrapolated luminosity function of 
quasars introduced in the simulation is shown in 
Figure 1 for z = 10; the finite box size causes 
the discreteness effects observed at high luminos- 
ity end. Our simulation box includes quasars up 
to a luminosity of about l0 12 Ls un at z = 10 (al- 
though quasars this luminous are rare, and there- 
fore appear in only a few lines of sight). 

In our simulation we include biasing of sources 
of ionizing and Ly-a radiation. Specifically, we 
assume that galaxies are biased with a bias factor 
of two, while quasars are biased with a bias factor 
of three. In this work we do not consider rcdshift- 
dependent bias factors. 

The specific simulation used in this paper has 



2 



10-' 



10- 6 - 
10- 7 - 




10 7 10 s 10 9 10 10 10" 10 12 

L [L su> ] 

Fig. 1. — Our adopted quasar luminosity function 
at z — 10, as extrapolated from Schirber and Bul- 
lock (2003). The dash-dotted line shows the lumi- 
nosity of a typical quasar which is expected to be 
present in almost all 10' observational beams. 

128 3 cells with an average comoving cell size of 
2/i _1 Mpc. The box size of 256/i~ 1 Mpc allows us 
to resolve individual H II regions around typical 
quasars at z ~ 8 — 10, although it is still too 
small to include the brightest high-redshift quasars 
found in the Sloan Digital Sky Survey (Fan et 
al. 2003). But since we are interested in typical 
(rather than the rarest, most extreme) features in 
the redshiftcd 21cm signal, this simulation is suf- 
ficient for our purpose. 

3. 21 cm "Spectral Dips" due to Quasar 
H II Regions 

We now investigate the characteristics of H II 
regions as they appear in the rcdshifted 21cm 
spectra from neutral hydrogen. The intergalactic 
hydrogen can be observed in absorption or emis- 
sion against the CMB, depending on whether the 
spin temperature is lower or higher than the CMB 
temperature. Although an absorption signal is ex- 
pected when the earliest emission of Ly-a photons 
occurs (which couples the spin and kinetic tem- 
perature of hydrogen, bringing the spin tempera- 
ture below the CMB temperature; see Madau et al. 
1997), the 21 cm signal in our simulation quickly 
turns into emission, as X-rays, far-ultraviolet ra- 
diation, and weak shocks from structure forma- 
tion heat the atomic gas to a temperature higher 
than the CMB (e.g., Chen & Miralda-Escude 2004; 



Gnedin & Shaver 2004). For the redshift range 
that we consider in this paper (corresponding 
to frequencies largely above the FM bands), the 
atomic gas is heated above the CMB tempera- 
ture and produces an emission signal. Therefore, 
H II regions should produce dips in the spectrum 
because they contain less neutral hydrogen than 
their environment, so their emission is reduced 
compared to the regions around them. We shall 
generally refer to these H II region features as dips. 
Note that these dips look like absorption features 
even though they have nothing to do with absorp- 
tion. 

The basic observable properties of the spectral 
dip of an H II region are its width and depth. The 
depth depends on the ionization fraction, gas den- 
sity, spin temperature, and also on the angular 
size of the H II region if this angular size is not 
resolved by the beam of the radio telescope that is 
used (the angular size of a typical observable H II 
region is about 1.5 — 2arcmin). The width of the 
spectral dip in frequency space is determined by 
the radial extent of the ionized region (assuming 
that peculiar velocities are small compared to the 
Hubble expansion across the H II region). 

In practice, the detected cosmological H II re- 
gions will probably be smaller than the angular 
size of the radio beam. However, if the H II re- 
gions are resolved in frequency, their frequency 
width determines their size. Thus, the average lin- 
ear (and angular) size for a sufficiently large sam- 
ple of observed H II regions can be estimated (for 
a given cosmological model) from the frequency 
data alone, and, taking into account the dilution 
by the known radio beamsize, the actual average 
depth of the spectral dip can be inferred from the 
measured dip depth. 

In this paper, however, we do not investigate 
such a measurement any further, because it will 
most likely require a large number of independent 
lines of sight. 

The average depth of a spectral dip is related 
to the difference in the amount of neutral hydro- 
gen present (on average) inside and outside an 
H II region. The observational determination of 
this quantity would introduce a new tool for test- 
ing models of reionization. In the H II regions 
the fraction of gas that is still atomic depends on 
how much of the gas in collapsed, low-mass ha- 
los remains self-shielded, and how much has been 
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ionized and evaporated from halos (e.g., Shapiro, 
Iliev, & Raga 2004), a process that is important 
to determine the clumping factors of ionized gas. 
Outside the identified H II regions, the fraction 
of gas already ionized depends on the presence of 
hard photons and low-luminosity sources produc- 
ing unresolved ionized regions. Determining the 
dip depth distribution will require a large number 
of independent lines of sight, and the theoretical 
interpretation will be complex. In the rest of this 
paper we focus on how spectral dips due to H II 
regions can be identified. 

In the absence of H II regions, the distribution 
of brightness temperature of the redshifted 21cm 
signal should be Gaussian, reflecting the imprints 
of cosmic density fluctuations, which are linear on 
the large scales that are most easily observable. 
The H II regions add additional features in the 
spectrum that are mainly caused by variations in 
neutral fraction. Gaussian regions of lower emis- 
sion correspond to regions of lower density, but 
H II regions usually have higher than average den- 
sity and a dramatically lowered neutral fraction. 
This can be used for distinguishing Gaussian fluc- 
tuations from those caused by the ionization from 
quasars in the simulation, but in real observations 
we will have to rely on the intrinsic properties of 
spectral dips to separate H II regions from density 
fluctuations. In the next two sections we will use 
our simulation to predict how such a separation 
can be made from the observational data. 

4. Analysis of Simulated Lines of Sight 

We now discuss how we use the output at con- 
secutive timesteps from the simulation described 
in §2 to create synthetic spectra of redshifted 
21 cm emission. These spectra reflect the large- 
scale distribution of gas density, spin temperature, 
and hydrogen neutral fraction. The spin temper- 
ature is computed self-consistently in the simu- 
lation including the effects of Ly-a pumping and 
collisions with electrons and neutral atoms (see 
Gnedin & Shaver 2004 for details). The pixel size 
in the spectra is set by the cell size of 2h~ 1 Mpc, 
which corresponds to a frequency range of about 
0.2 MHz at z ~ 9. 

The synthetic spectra are created by casting a 
random line of sight through the simulation box. 
Our simulated lines of sight span a significant frac- 



tion of the horizon distance, so we simulate the 
line of sight on a fixed past light-cone instead of 
a fixed cosmic time, using consecutive output files 
at different timesteps of the simulation. We also 
account for a finite observing beamsize by convolv- 
ing the 21cm signal over a plane perpendicular to 
the line of sight with a Gaussian beam, with a 
physical radius that varies over the line of sight 
according to a fixed angular size of a radio beam 

e 6 . 

The starting redshift of our spectra is z = 13, 
which corresponds to a frequency of 101 MHz. 
This redshift was chosen to place the resulting fre- 
quency in the range of next-generation radio ob- 
servatories. We used 450 independent lines of sight 
to accumulate sufficient statistics. 

After the synthetic spectra were created in this 
way, we analyzed the brightness temperature dis- 
tribution for signatures of H II regions. The pro- 
cess we followed for this analysis consists of the 
following steps. First, we fit the mean signal (Tb) 
with a fourth-degree polynomial in order to re- 
move variations on a very wide frequency scale. 
This process, similar to the continuum fitting used 
in the analysis of Ly-a forest spectra, is intended 
as an approximation to the removal of the com- 
bined galactic and extragalactic foregrounds and 
the cosmological mean signal in an observation. 
Next, we calculate the fluctuating signal as: 

AT B =T B -(T B ), (1) 

to obtain a fluctuating signal with zero mean that 
can be searched for H II region dips. For the analy- 
sis we disregard all parts of the spectrum in which 
the emission exceeds the mean signal, which corre- 
spond to regions of high gas density that are part 
of linear Gaussian fluctuations. We are then left 
with a spectrum showing dips due to both Gaus- 
sian fluctuations (corresponding to regions of low 
gas density) and H II regions around bright sources 
(corresponding to regions of very low neutral frac- 
tion compared to the average medium) . 

5. Results 

In Figure 2 we show an example of a typical line 
of sight with the mean signal subtracted. The top 
panel shows the brightness temperature fluctua- 
tion AT B versus v after subtracting the mean sig- 
nal. A significant fluctuation due to an H II region 
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Fig. 2. — Top panel: the fluctuation in the bright- 
ness temperature ATb as a function of frequency 
in a 10' beam with a 0.5 MHz bandwidth. Also 
shown are the 5cr sensitivity limits for a 1000 hrs 
integration time (Shaver et al. 1999). Bottom 
panel: the distribution of gas density (gray) and 
hydrogen neutral fraction (black) along the same 
line of sight. The H II region spectral feature can 
be seen at v = 144 MHz. 

at v = 144 MHz is easily identifiable in the spec- 
trum. Superimposed are the 5a sensitivity limits 
for an integration time of 1000 hrs and a band- 
width of 0.5 MHz (we discuss our choice for the 
bandwidth below). 

It is important to note here that our results dif- 
fer to a small degree from the predictions of Tozzi 
et al. (2000) and Wyithe & Loeb (2004) in that we 
do not observe an increase in the redshifted 21 cm 
emission just outside the H II region. This is due 
to the fact that the main source of Ly-a radia- 
tion that couples the gas kinetic temperature and 
the spin temperature of the 21 cm transition in our 
simulation is normal galaxies. Since, in our simula- 
tion, galaxies typically form before the quasars, an 
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Fig. 3. The auto-correlation function of the 
brightness temperature ATb as a function of 
FWHM frequency width 8v. 

H II region around a typical high-rcdshift quasar 
expands in the already pre-heated gas, in which 
the spin temperature of the 21 cm transition is al- 
ready much higher than the CMB temperature, 
so the extra Ly-a radiation from the quasar does 
not lead to an increase in the emission. We would 
like to re-emphasize here that, since we are treat- 
ing the radiative transfer in full 3D, we do include 
the proximity effect of the increased Ly-a radia- 
tion around quasars, but this effect, as our results 
indicate, is not significant. 

To distinguish H II regions from under-dense re- 
gions with low neutral hydrogen density, we now 
look at the distribution of density and neutral frac- 
tion in the same frequency range in the data ob- 
tained from the simulation. 

The bottom panel in Figure 2 depicts the vari- 
ations in these two parameters. The neutral frac- 
tion shows a dip at v = 144 MHz, whereas the 
density shows a definite increase. Consequently 
we can determine that this fluctuation results from 
a high density region containing more ionized gas 
than the average. We expect to find quasars in 
high density regions since their population is bi- 
ased. If the radiation intensity were approximately 
uniform, we would expect to find lower densities 
of neutral gas in under-dense regions of the IGM, 
because of the decreased recombination rate. Here 
the high density and ionized fraction of the IGM 
allow us to conclude that the spectral dip is caused 
by an H II region surrounding a bright quasar. 

The fluctuations in the brightness temperature 
on large scales are expected to be Gaussian, re- 
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Fig. 4. — The histogram of ATb for the normal- 
ized spectrum from Figure 2: notice the long tail of 
the non-Gaussian contribution from H II regions. 

fleeting the large-scale variations in the cosmic gas 
density and neutral fraction. However, on suffi- 
ciently small scales the fluctuations become corre- 
lated. Thus, unless the beam size of a radio tele- 
scope is small enough (less than about 1 arcmin) 
to resolve individual H II regions around high red- 
shift galaxies, one would expect the cosmologi- 
cal signal to become smooth on scales comparable 
to the correlation length of high redshift galax- 
ies. Figure 3 shows the auto-correlation length of 
the brightness temperature from the line of sight 
shown in Fig. 2 Thus, decreasing the bandwidth of 
observation to below about 0.5 MHz will not result 
in measuring new information in the signal unless 
it also coincides with a decrease in the beam size 
(this is also apparent from Fig. 11 of Gncdin & 
Shaver 2004). In the rest of this paper, we assume 
a fixed value of 0.5 MHz for the bandwidth for our 
fiducial beam size of 10'. 

In addition to investigating individual spec- 
tra as shown above, we can use the variations 
of brightness temperature around the average to 
look at the distribution of spectral dips statisti- 
cally. Figure 4 shows a histogram of the difference 
in brightness temperature of the spectrum in Fig- 
ure 2. One can see the tail of the distribution 
on the negative ATg side that correspond to dips 
in the spectrum due to quasar H II regions. The 
peak centered around ATb = is mainly due to 
the Gaussian fluctuations, but the long tail on the 
negative ATb reaching ATb = —4.5 mK shows the 
significant contribution from regions that are due 
to ionized regions and not from Gaussian density 
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Fig. 5. — Top panel: in black is shown the distri- 
bution of equivalent widths and sizes of the spec- 
tral dips from 450 lines of sight. The gray data 
shows the same distribution for a pure Gaussian 
distribution. Bottom panel: A similar plot but for 
equivalent width versus the depth of the spectral 
dip distribution. 

fluctuations. 

We use a simple calculation of the equivalent 
width of each dip as a way to measure the dip 
depth. We integrate the area below the mean 
signal for each dip (i.e. a connected region with 
ATb < 0), and compile a distribution for all the 
dips of 450 different spectra. This distribution is 
shown in Figure 5 by plotting the equivalent width 
versus the dip width in frequency space (measured 
as FWHM). As a comparison, we also show the 
same distribution for artificial lines of sight that 
contain no H II regions, so that any dips are en- 
tirely due to Gaussian variations in cosmic density. 

It is clear that the largest equivalent widths and 
depths correspond to H II regions around quasars, 
since the Gaussian density variations do not gen- 
erate such large dips. The radial sizes, as shown 
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by Av>, of both samples are comparable, since the 
Gaussian density fluctuations also yield the wide 
dips seen in the spectra. However, the Gaussian 
fluctuations are much shallower, as can be seen by 
the smaller values in the equivalent width. There 
is also a range of smaller equivalent widths that 
is not reached with purely Gaussian fluctuations. 
This could be due to faint quasars in the simula- 
tion causing the small features. Judging from this 
graph, all spectral dips with equivalent width in 
excess of 3mKMHz or depth in excess of 1.5 mK 
are caused by individual H II regions. 

We can also look at cumulative distributions 
of various parameters characterizing these spec- 
tral dips, which are shown in Figures 6, to de- 
termine which are due to Gaussian fluctuations 
compared to the ones due to H II regions. The fig- 
ure shows distributions for the simulated fluctua- 
tions in black and contrasts it with the distribution 
for pure Gaussian fluctuations in gray. From Fig- 
ure 6a it is again clear that all spectral dips with 
equivalent widths in excess of about 3 mK MHz or 
depths in excess of about 1.5 mK are coming from 
high-rcdshift H II regions. The widths of the dips 
in frequency space, however, cannot be used to 
separate H II regions from linear Gaussian fluctu- 
ations. From these graphs we can infer that there 
is about a one in three chance to detect an H II re- 
gion with equivalent width in excess of 3 mK MHz 
in each line of sight, but narrower deep dips (with 
depth in excess of 1.5 mK) are present in almost 
any line of sight. 

The results shown in Figures 6 were obtained 
for a bandwidth of 0.5 MHz, but as long as the 
bandwidth is less than about 1 MHz, spectral dips 
from H II regions are resolved (as can be seen from 
Fig. 5a), and the distributions of equivalent widths 
and depths remain the same. 

6. Discussion 

We have shown that by far the easiest to ob- 
serve features of the cosmological signal of red- 
shifted 21cm emission from neutral hydrogen 
in the pre-reionization era are the spectral dips 
due to individual H II regions from high-redshifts 
quasars. The largest of these features can be eas- 
ily distinguished from Gaussian fluctuations due 
to large-scale structure: for example, we find that 
all spectral dips with equivalent width in excess of 
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Fig. 6. — Top: The cumulative distribution of 
spectral dips per line of sight as the function of 
the dip width Az/ (measured as FWFM - top), the 
equivalent width (middle), and the depth (bottom) 
for the simulated spectra (black) and for purely 
Gaussian fluctuations (gray). 
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3 mK MHz or depth in excess of 1.5 mK come ex- 
clusively from quasar H II regions (for a 0.5 MHz 
bandwidth) . 

These dips are not rare - for our adopted plausi- 
ble (albeit uncertain) extrapolation of the quasar 
luminosity function, detectable dips are present 
in almost every line of sight. This conclusion 
may sound counterintuitive to the common wis- 
dom that "high redshift quasars are rare", but 
we need to emphasize that the quasars that cre- 
ate detectable H II regions are much less luminous 
(and much more numerous) than, say, SLOAN 
quasars. Besides, a 10' beam (for example) go- 
ing from z = 6 to z = 12 (the onset of the FM 
radio band) contains a significant volume of about 
500,000(/i _1 Mpc) 3 . According to our model of 
the quasar luminosity function, there is about one 
quasar within this volume with blue rest-frame lu- 
minosity of 2 x l0 10 Lsun, and this luminosity is 
sufficient for a typical quasar to create an H II re- 
gion that can be distinguished from Gaussian fluc- 
tuations due to large-scale structure. For a differ- 
ent model of the quasar luminosity function, our 
results can be simply scaled with the abundance 
of quasars above this characteristic luminosity. 

Knowing the shapes and sizes of the ionized re- 
gions around quasars at high redshift could help 
to determine better the physical processes near 
super-massive black holes at the very first stages of 
structure formation in the universe. In particular, 
if observed spectral dips in the redshifted 21cm 
emission can be identified with specific quasars ob- 
served in the infrared by, say, JWST, then a sig- 
nificant constraint can be placed on the lifetimes 
of quasars as a function of their luminosity. Such 
a constraint is difficult to obtain using other meth- 
ods. 

The depth distribution of spectral dips caused 
by H II regions can put a constraint on the resid- 
ual neutral fraction inside H II regions. This neu- 
tral fraction would likely be locked in high-density, 
self-shielded gas clumps in collapsed halos (ob- 
served as Lyman limit systems in ultraviolet spec- 
tra), which play the dominant role for determining 
the clumping factor of ionized gas during rcioniza- 
tion. 

Under the assumption that the distribution of 
H II regions is statistically isotropic, it will be pos- 
sible to put constraints on the residual neutral 
fraction inside H II region, and, since this neutral 



gas most likely is locked inside Lyman-limit sys- 
tems, an important constraint can be put on the 
abundance of low mass objects. 

In conclusion, the spectral features predicted 
in this paper - the spectral dips due to H II re- 
gions around individual high-rcdshifts quasars - 
will be the most prominent signals in the red- 
shifted 21cm line of neutral hydrogen from the 
pre-reionization era. Their spectral sharpness will 
easily distinguish them from the smoothly vary- 
ing continuum spectra of the foreground emissions, 
and their strength should make them detectable in 
reasonable integration times with the new genera- 
tion of low frequency radio telescopes. 
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